Unit Wspli;

interface

uses Objects, WinTypes, WinProcs, WinDos, Wincrt, BWCC, Strings,
     OStdDlgs, OWindows, ODialogs, OMemory,
     CartoCmds, tit_win;

CONST maxr = 256;
      maxc =  3;

TYPE
    op_type = (plus,minus,multiply,divide,invert,cmod);
    short_vector  =          ARRAY[0..maxc] of real;

    plong_vector  =          ^long_vector;
    long_vector   =          ARRAY[0..maxr] OF real;

    precmat       =          ^recmat;
    recmat        =          ARRAY[0..maxr,0..maxc] OF real;

function cubic_splines(x,y: plong_vector; n: integer; var s: recmat):integer;
procedure calc_spline(xv: long_vector; coef: recmat;
                      n: integer; x: real; var y: real);

implementation

function cubic_splines(x,y: plong_vector; n: integer; var s: recmat):integer;

const  nearzero = 1.0e-20;
var i: integer;
    eqn: precmat;
    bb,cc,hh,ff,fi,hi,nextf: real;

    h,bvec: plong_vector;
label 999;

procedure tridiag(mat: recmat; n: integer; var vec: long_vector);
const  nearzero = 1.0e-20;
var i: integer;
    pivot,mult,ci,bi: real;
begin
  for i := 1 to n-1 do
  begin
    pivot := mat[i,3];
    ci := mat[i,2];
    bi := mat[i,1];
    mult := mat[i+1,0] /pivot;
    if abs(mult) > nearzero then
    begin
      mat[i+1,0] := mult;
      mat[i+1,3] := mat[i+1,3] - mult * ci;
      mat[i+1,1] := mat[i+1,1] - mult * bi;
    end
    else mat[i+1,0] := 0;
  end;
  vec[n] := mat[n,1] / mat[n,3];
  for i := n-1 downto 1 do
    vec[i] := (mat[i,1] - mat[i,2] * vec[i+1]) / mat[i,3];
end;

begin
   cubic_splines := -1;

   h := nil; bvec := nil; eqn := nil;
   ier := _getmem('-',pointer(h),sizeOf(long_vector));
   if ier <> 0 then goto 999;
   ier := _getmem('-',pointer(bvec),sizeOf(long_vector));
   if ier <> 0 then goto 999;
   ier := _getmem('-',pointer(eqn),sizeOf(recmat));
   if ier <> 0 then goto 999;

   for i := 0 to n-1 do h^[i] := x^[i+1] - x^[i];
   hh := h^[0];
   ff := y^[0];
   nextf := y^[1];
   for i := 1 to n-1 do begin
     fi := nextf;
     nextf := y^[i+1];
     hi := h^[i];
     eqn^[i,1] := 3 * ((nextf - fi)/hi- (fi - ff)/hh);
     eqn^[i,3] := 2 * (hh + hi);
     eqn^[i,0] := hh;
     eqn^[i,2] := hi;
     hh := hi;
     ff := fi;
   end;
   tridiag(eqn^,n-1,bvec^);
   bb := 0;
   hh := h^[0];
   s[0,2] := bb;
   s[0,0] := y^[0];
   cc := (y^[1] - s[0,0])/hh - bvec^[1] * hh/3;
   s[0,1] := cc;
   for i := 1 to n - 1 do
   begin
     s[i,2] := bvec^[i];
     s[i,0] := y^[i];
     cc := (s[i,2] + bb) * hh + cc;
     s[i,1] := cc;
     s[i-1,3] := (s[i,2] - bb)/(3 * hh);
     bb := s[i,2];
     hh := h^[i];
   end;
   s[n-1,3] := -bb/(3 * hh);

   cubic_splines := 0;
999:
   _freemem('-',pointer(h), sizeof(long_vector));
   _freemem('-',pointer(bvec), sizeof(long_vector));
   _freemem('-',pointer(eqn), sizeof(recmat));
end;

procedure calc_spline(xv: long_vector; coef: recmat;
                      n: integer; x: real; var y: real);
var i,j,k: integer;
    coef_vector: short_vector;
    dif: real;

function poly_calc(x_in: real; coef_vector:short_vector;
                    order: integer): real;
var i: integer;
    y_out : real;
begin
  y_out := coef_vector[order];
  for i:= order - 1 downto 0 do
    y_out := y_out * x_in + coef_vector[i];
  poly_calc := y_out;
end;

begin
  i := -1;
  repeat
    i := i + 1;
  until ((x >= xv[i]) and (x < xv[i+1])) or (i = n);
  dif := x - xv[i];
  for j := 0 to 3 do
    coef_vector[j] := coef[i,j];
  y := poly_calc(dif,coef_vector,3);
end;

end.

